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Abstract 

A domain integral method employing a specific Green's function (i.e., incorporating some 
features of the global problem of wave propagation in an inhomogeneous medium) is developed 
for solving direct and inverse scattering problems relative to slab-like macroscopically inhomo- 
geneous porous obstacles. It is shown how to numerically solve such problems, involving both 
spatially-varying density and compressibility, by means of an iterative scheme initialized with 
a Born approximation. A numerical solution is obtained for a canonical problem involving a 
two-layer slab. 

1 Introduction 

This work was initially motivated by two problems: i) the design problem connected with the 
determination of the optimal profile of a continuous and/or discontinuous spatial distribution of 
the material/geometric properties of porous materials for the absorption of sound [Zj and ii) the 
retrieval of the spatially-varying mechanical and geometrical parameters of bone for the diagnosis 
of diseases such as osteoporosis [Tflj . 

Such inverse problems [2U] can be decomposed into two sub-problems: i) the determination of the 
constitutive and conservation relations linking the various spatially-variable mechanical parameters 
of the porous medium to its response to an acoustic solicitation, and ii) the resolution of the 
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wave equation in an inhomogeneous porous medium (for instance, within the Biot, or rigid frame 
approximations). Here we focus on the second point. 

In ^2]) it is shown that the wave equation describing the propagation in a macroscopically- 
inhomogeneous porous medium in the rigid frame approximation can formally take the form of 
the usual acoustic wave equation in a macroscopically-inhomogeneous fluid (in which the micro- 
scopic features of the porous medium are homogenized) with spatial (and frequency) dependent 
compressibility /c e (x, uS) and density p e (x,u). 

The present work deals with a method of resolution of direct problems involving acoustic wave 
propagation in a macroscopically-inhomogeneous fluid medium, whose density and compressibility 
are both space dependent, this being a prerequisite to the resolution of related inverse problems. 

This topic is also of great interest in quantum physics (inverse potential scattering |2ll37 |l38| l39|). 
ocean acoustics El EU EEU HUE EES] (detection of inhomogeneities, sediment exploration, influence 
of seawater and seafloor composition and heterogenity on the long-range propagation of acoustic 
waves in the sea, ...), seismology [Hl4()|H2*] (determination of the internal structure and composition 
of the Earth via seismic waves,...), geophysics |421 1251 l4"51 144j (characterization of soil, detection of 
geological features such as hydrocarbon reservoirs, ...), optics and electromagnetism [471 141 j (design 
and characterization of materials having specified response to waves, detection of flaws,...). 

The wave equation in an inhomogeneous medium can be solved in a variety of manners: via the 
wave splitting method [313 El E01 > the transfer matrix method HFJIH (for piecewise constant media), 
integral methods |2D | I39 | 137] . or purely numerical (e.g., finite-element |19| or finite-difference jS]) 
methods. The methods dedicated to inverse problems are wave splitting and linearisation |381 139j 
techniques deriving from the integral formalism. The two most widely-known approximations for 
the Fredholm equations of the second kind involved in the integral formalism (at least when the 
density is constant in the acoustic context) are the Born approximation |35l I39j and the Rytov 
approximation [23! ■ We will focus on the Born approximation, despite the fact that several authors 
|27l I13j have shown that the Rytov approximation is valid under a less restrictive set of conditions 
than the Born approximation. 

We postulate, and show, that the accuracy of the Born approximation can be increased by 
the use of the integral formulation together with a specific Green's function |37 | I39[ H6] . In most 
of the articles dealing with the Born approximation and other linearisation methods, the prob- 
lems are often simplified by considering the density to be constant. Herein, we consider both the 
compressibility and the density to be spatially- variable. This induces supplementary difficulties, 
because it can lead to meaningless integrals (involving first and/or second space derivatives of the 
density), especially when the variation of the density is not continuous, and/or because it requires 
the evaluation of the first space derivative of the pressure field. We will show how to deal with 
these problems. 

The usual first-order Born approximation is an outcome of the integral formulation employing 
the free-space Green's function (FSGF) and consists in approximating the pressure field in the 
integrand (corresponding to the pressure field inside the heterogeneity) by the field in the absence 
of the heterogeneity (i.e. the incident field), this being equivalent to the asssumption that the 
diffracted field is negligible compared to the incident field. Although this method usually provides 
good results for small contrasts between the mechanical parameters of the inhomogeneity and those 
of the host medium, its accuracy decreases in the case of dissipative media and larger contrasts. 

Moreover, iterative schemes initialized with the zeroth-order Born or Rytov approximations 
often diverge in practice. This difficulty can be partially resolved by employing the modified or 
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distorted Born approximation j^H] which basically consist in acting on one of the three terms 
involved in the integrand (the Green's function, the contrast function and the field inside the 
heterogeneity) . 

The method described herafter allows us to act on all three terms. The central idea of our 
method consists in reducing the eigenvalues of the kernel of the integrand (thought to be the 
cause of the difficulties with the usual iterative Born scheme) by employing a Green's function-the 
so-called Specific Green's function (SGF)— of a canonical problem which is close, in some sense, 
to the original problem. The specification of the initial solution was already treated in [1^1 in 
connection with the resolution of an inverse problem. The chosen problem (also canonical) is that 
of the diffraction of an incident plane wave, propagating in the host medium, by a two-component 
slab (each component being a homogeneous layer) considered as a single inhomogeneous slab. In 
the present instance, the close canonical problem involves a slab filled with a macroscopically- 
homogeneous fluid-saturated porous medium surrounded by the same fluid (air) medium as the 
original macroscopically-inhomogeneous fluid-saturated porous medium. 

We will show, for this example: i) how to construct an appropriate specific Green's function 
(SGF), ii) how to incorporate the latter into the integral formulation, and iii) how the resulting 
integral equation can be solved by an iterative scheme initialized by a modified Born approximation. 

The results are compared to those of the analytic solutions (obtained by the transfer matrix 
method (TMM)) for the two-component slab and found to be in good agreement with the latter, 
both in transmission and reflection and for several angles of incidence. The iterative scheme initial- 
ized with the modified Born approximation converges rather rapidly, i.e., within 5 to 7 iterations for 
our example. This demonstrates the efficiency of our SGF interative scheme for the resolution of the 
direct problem relative to wave propagation in the presence of a macroscopically-inhomogeneous 
fluid-saturated porous slab. 

2 Use of the specific Green's functions in the domain integral 
formulation to solve direct scattering problems 

2.1 An example of a direct scattering problem 

The type of direct problem we deal with is illustrated in figure^,. This problem involves spatially- 
dependent compressibility and density. As will be shown further on, the spatial variability, and 
discontinuity, of both these quantities, and in particular of the latter, can produce some difficulties 
in the domain integral formulations (but not in the TMM formulation). The wave equation for 
such problems is given in appendix ^ to solve them in optimal manner, we treat the auxiliary 
problem depicted in figure^. 
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Figure 1: Configuration of the direct problem of a.: a fluid heterogeneity within a fluid- like slab, 
b.: diffraction of a wave radiated by applied sources from a fluid slab of boundaries T a and r& 
immersed in a fluid. 

2.2 Specific Green's function corresponding to the propagation of waves radi- 
ated by interior and exterior line sources in the presenece of a homogeneous 
fluid-like layer immersed in a homogeneous fluid host medium 

2.2.1 Features of the problem 

The sagittal plane (cross-section) view of the scattering configuration is given in figure \T]p. As we 
are dealing with a Green's function, the supports of the sources reduce to dots in the figure, i.e., 
the sources are line sources. The homogeneous fluid-like layer is oriented horizontally (i.e. the 
normal to both of its faces is along the X2 axis); its thickness is I, and the medium M 1 therein is 
homogeneous. The geometry and composition of the layer are thus invariant with respect to x%. Q,\ 
designates the trace of the layer in the x\ — X2 cross-section plane. T a and r& designate the traces 
of the lower and upper faces respectively of the layer in the x\ — X2 cross-section plane. The unit 
vectors normal to r a and are designated indistinctly by v. The X2 coordinates of r a and are 
designated by a and b respectively. 

The layer is immersed in a (host) fluid M°. The trace of the host medium domain below (above) 
the layer in the x\ — X2 cross-section plane is designated by Q + (Sl - ) • 

The (direct scattering) problem is to determine the response g 0+ within f?o+> 9° within f2 - 
and g 1 within Qi for line sources that are located either within Oo+> ^o- or This response 
constitutes the specific Green's function we are looking for. 

Let y designate the vector from O to the location of the line source. The Green's function in 
Vtj is designated by g J (x,y), which means the response at x due to line sources located at y. 

2.2.2 Governing equations 

Rather than to solve directly for g 3 (x,y,t), we prefer to deal with its Fourier transform g J (x,y,uj) 
defined by: 

/oo 
^(x,y, W )e- 1 ^;i = 0+,l,0- (1) 
-co 
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The mathematical translation of the boundary- value problem in the space- frequency domain is: 

A+(k j ) 2 ]g j (x,y,u;) =-«5(x-y) : Vx G j = 0+ 1, 0", y G O 0+ , ^ or tt - , (2) 



Vx G r a { 



Vx G r b < 



5° (x,y,w) -5 1 (x,y,w) = 

-oi/(x) • V 5 0+ (x, y, w) - — i/(x) • V*? 1 (x, y, w) = 
5 1 (x,y,w) - g° (x,y,w) = 

-^i/(x) • V3 1 (x,y,w) - ^iv(x) • V<7°~ (x,y,w) = 



(3) 



(4) 



g- 7 (x, y,u>) — G- 7 (x, y, u) ~ outgoing waves, ; Vx G flj ; j = + , 1,0 , ||x|| — > oo (5) 
wherein G 7 is the free-space Green's function in the medium M- 7 given by 

G»(x,y,u;) = ^Vdl* - y||) = A |°° exp (i^ (x x - yi ) + \k{ \x 2 - y 2 \) ^ , (6) 



with i?Q 1; the zeroth-order Hankel function of the first kind, k? 2 = \J (A;- 7 ) 2 — {k\) 2 such that 
ft fe) > and 9 fe) > 0, j = 0, 1 for to > 0. 



2.2.3 Field representations 

We shall henceforth: i) drop the cj-dependence with the understanding that it is implicit in all the 
field functions and ii) employ the cartesian coordinates (xi,x 2 ) of x and (2/1,1/2) of y. 
We use the separation of variables technique to obtain 

/OO Jl, 
E 0+ exp(ita + ifc 2 °(z 2 -a))-/ (7) 
-00 fc 2 

(A 1 exp (-ik\x 2 ) + B 1 exp (ik\x 2 )) exp (ik lXl ) -j- (8) 
-00 ^2 

A°"exp(iA;ixi-i^(ar 2 -6))-/ (9) 

-00 &2 

wherein Hq . is the Heaviside function 



Ho 7 (y) 



1 if y G % 

if y G fii, i / j 



(10) 
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2.2.4 Application of the transmission conditions 

In cartesian coordinates and on account of the orientation of the two faces of the layer: 

i/(x)-V^ = ^-^. 

ox 2 



(11) 



After introducing the fields expressions, eqs.0, (jHj) and Q into the boundary conditions eqs.© 
and (jlj), we multiply these relations by exp (— \K\X\) and then integrate from —00 to +00, using 
the identity 



f 

j — < 



exp (i {k\ — K\) x\) dx\ = 2ir5 {k\ — K\) , 5 (k\ — K\) being the Kronecker symbol (12) 



so as to obtain the matrix equation 
/ 



1 
1 





_ k O e -ikla 


_ A ,0gifcia 


k\ 


k\ 


pOg-ifc^a 




P 1 


P 1 


_ fc 0g-lfci& 




k\ 


k 2 


_ ( p e -\k\b 




P 1 


P 1 






1 
1 



A 1 
B l 

A " 



le 



i*S(l/2-a)i 



47T 



g^ 2Vy ,- u ; HQ 



Ml p l 



_ e i^-y 2 ) HfV + ^gi^fc-^H^ 



Hn _ + V^ 2 ~ 6)H ^ 



(13) 



2.2.5 Final expressions of the specific Green's function 

Once the matrix system (|13[) is solved for B 0+ , A 1 , B 1 and A , and these expressions are introduced 
into the expressions of the fields (J7J), (jSJ), @, we get: 

5 0+ (x,y) 



1 JiMzi-ai)+ifcg|2 2 -2/2 n " dfc 



47T 



,[ifci(xi-j/i)+ifc^ 2 ] 



, 4vr (2a°a 1 cos {k\l) - i ((q ) 2 + (a 1 ) 2 ) sin 
>°(y2-2a) gin ^ ((a } 2 _ (q 1 )2 ) lgo+ + 2ie -ifc 2 °(^+0 a + 

2ie~ ife ° V (a 1 cos (y 2 - b)) - ia° sin (** (y 2 - 6))) ^ 

K 2 



dh , (14) 
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dk-\ 



iki(xi—yi) 



4vr (2a°a 1 cos (fc^) - i ((a ) 2 + (a 1 ) 2 ) sin 



Hr 



A, 2 



2e i*S(io-a) ^ a i a o CQS ^1 (x2 _ fo)) _ i ( a 0) 2 sin ( X2 _ 



-+ 



-+ 



a 1 ) 2 — (a ) 2 ) cos {k\ (^2 + 2/2 — a — 6)) + exp (ik\l) (a — a 1 ) 2 cos (k\ (x2 — y2 



He 



A<2 



(15) 



,[ifci(a;i-j/i)-ifc^X2] 



3 4vr (2a°a 1 cos (jfe£l) - i ((a ) 2 + (a 1 ) 2 ) sin 

2 ie -i*S(»-io) a i a o5^L + e -ife 2 (?/2-2 b ) gin ^i^ ^ (a o )2 _ (a i )2 ) 

fc 2 fc 2 



2ie lfe 2 6 a 1 (a 1 cos (a - y 2 )) 



ia° sin 



dfci . (16) 



2.3 Use of a specific Green's function to solve the direct problem of pressure 
wave scattered by an inhomogeneous fluid-filled slab 

We treat the 2D fluid acoustic direct problem illustrated in figure El I n the absence of the het- 
erogeneity, occupying the domain O2, the configuration is that of the closed layer domain S7i 
occupied by a known homogeneous fluid M 1 with (spatially-constant) acoustic parameters (k , p 1 ), 
surrounded by the open domain Qo occupied by a known homogeneous fluid M° with (spatially- 
constant) acoustic parameters (fc°, p°). 

In the presence of the heterogeneity, localized to the domain 0,2 £ the problem is to solve the 
scattering problem for spatially-varying acoustic parameter functions (fc 2 (x), p 2 (x)) of the medium 
M 2 filling in the subdomains Qq+ an d ^0- when the slab is probed by an incident wave. 



7 




Figure 2: Configuration of the direct problem of a fluid heterogeneity within a fluid-like slab. 



2.3.1 Governing equations for scattering from a heterogeneous layer, included be- 
tween T a and ]?{,, probed by a cylindrical wave radiated by a cylindrical source 
whose support is 



Let f^i = r?i U Og. Then the governing equations for the pressure field are: 

[A + (fc ) 2 ] p 0+ (x)=- S (x);xGtt 0+ , 



A + (A;(x)) 2 l p\ X ) = • Vp^x); x G n is 



P(x) 



fc(x) 



P(x) 



A' 1 




; x G Cli 


A; 2 


(x) 


; x G fi 2 


P 1 




; x G Oi 


P 2 


(x) 


; x G S7 2 


P°~ 


(x) 


= 0; x G Q - 


P*( 


x) = 


0; x g r , 



(17) 
(18) 
(19) 

(20) 

(21) 
(22) 
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li/(x) • Vp 0+ (x) - li/(x) • V/(x) = 0; x G r o , (23) 

p 1 (x)-p°"(x) = 0;xGr fe , (24) 

-^(x) • V/(x) - li/(x) • Vp°» = 0; x G r 6 , (25) 

p 0+ (x), p*(x) and p° (x) ~ outgoing waves, ||x|| — > oo (26) 
The previously-given governing equations for the specific Green's function can be rewritten as: 

[A + (K(x)) 2 ] 0(x, y) = -<5(x - y); x 6 I 2 , y e R 2 (27) 

f k° ; xeH 0+ U Q - 

K(x) = . (28) 
{ k ; x G Hi 

f <?°(x,y) ; x£fi 0+ U ft - 

S(x,y)= (29) 
I s (xj) ; x G fii 



,0H 



(x,y)-< 7 1 (x,y)=0;xGr a , (30) 



^i/(x) • V 5 0+ (x, y) - li/(x) • V^x, y) = 0; x G T a , (31) 

<7 1 (x,y)-< 7 °"(x,y)=0;xGr b , (32) 

-^(x) • Vg l (x, y) - li/(x) ■ V<7°~ (x, y) = 0; x G r 6 , (33) 

g(x, y) ~ outgoing waves, ||x|| — > oo (34) 

2.3.2 Towards a domain integral representation of the pressure field in Q + 

In obvious short-hand notation (in addition: d v := v ■ V), we obtain from the previous governing 
equations: 

^ [A + (Jb°) 2 ] = -/V ; in fl,* (35) 

p 0+ [A + (k ) 2 ] g 0+ = -p 0+ 5 ; in (l 0+ (36) 
so that integrating the difference of these two equations over Q + , we obtain 

I (g 0+ Ap 0+ -p 0+ Ag 0+ )dn = -[ g 0+ s°dn+[ p 0+ 5dQ (37) 
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or, after use of Green's theorem and the sifting property of the 5 distribution: 

/ (<? 0+ d„ P 0+ - P 0+ d^ + ) d 1+ 

/ (g 0+ d uP 0+ - P 0+ d u9 0+ ) d 7 + / g 0+ s dn = P 0+ (y)n n0+ ( y ) 

We develop, for the domain integral representation of this pressure field, the integration over rg+, 
so that to obtain: 



d 7 



[ (g° + d„p 0+ - p 0+ d„g 0+ ) dn= ( g 0+ \d u p Q+ - ifcy + j d 7 - / P° + \du9 0+ - ik°g 0+ 

Jpoo V / Jpoo L J Jpoo L 

0+ 0+ 0+ 

(39) 

It is readily shown that both of the integrals on the right hand side of this expression vanish due 
to the fact that both p 0+ and g 0+ satisfy the (frequency domain) rediation condition at infinity. 

2.3.3 Towards a domain integral representation of the pressure field in Q\ 

In obvious short-hand notation, we obtain from the previous governing equations: 

[A + (fc(x)) 2 ] p 1 = ^ • Vp 1 ; [A + (k 1 ) 2 ] p 1 = [(k 1 ) 2 - (k(x)) 2 } p 1 + ^ • Vp 1 = -<t(x) (40) 
Consequently, 

g 1 [A + (A: 1 ) 2 ] p 1 = -g l o5 ; in (41) 

p 1 [A + (A; 1 ) 2 ] g 1 = -p l 5 ; in n ± (42) 

so that, integrating the difference of these two equations over Q\, and after use of Green's theorem 
and the sifting property of the 5 distribution, we obtain: 

- / tfdvP 1 - p'd.g 1 ) d 7 + I (g l d v p x - p l d u g l ) d 1 + / g 1 od£l = p\y)^ (y) (43) 

JT a JT b JQ! 

which yields, on account of the transmission conditions: 

-/ (g 0+ d u p 0+ -p 0+ d u g 0+ )d 1 + f ( 5 o-^°--A^°-)d 7 + 4/ 9^dn = 

^p 1 (y)H„ 1 (y) (44) 

2.3.4 Towards a domain integral representation of the pressure field in Q - 

In obvious short-hand notation, we obtain from the previous governing equations: 

g°- [A + (k ) 2 } p°- = ; in n _ , (45) 
p°~ [A + (At ) 2 ] g°~ = -p°~ 5 ; in Q - , (46) 
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so that, integrating the difference of these two equations over O -, and following the procedure 
used in the last two subsections, we obtain: 



g°-d v p°- - p°~d u g°- ) d 7 = P u ~ (y)Hoo- (y) • 



(47) 



2.3.5 Domain integral representations, without boundary terms, of the pressure fields 
in O + , $7i and f2 - 

The addition of PS]), ijlljl and (gTJ) gives 

f g o + s ° d n + £ [ g l adn=p 0+ (y)R n0+ (y)+p\y)R nl (y)+p 0+ (y)R n0+ (y) , (48) 
Jn 0+ P 

from which it ensues, on account of the properties of the domain Heaviside function: 



P° + (y)=/ 9 0+ s dn + ^l g\x,y) (k(xf - (fc 1 ) 2 ) 



o+ 



P JQi 



p x (y) = 4 [ g 0+ s°dn+ [ gH^y) 

P Jn 0+ JUj 



Vp 



Vp 



• V 



(fc(x) 2 -(^) 2 )-^-V 
P 



p^x^O, Vy G fV » (49) 



p 1 (x)dn, Vy G ill , (50) 



P°"(y) 



g^s^n+Pj g\x,y) 
n 0+ P Jfi! 



(fe(x) 2 - (fc 1 ) 2 ) 



Vp 



V 



^(xjdfi, Vy G fi„- . (51) 



2.3.6 Other integral representations, without boundary terms, of the pressure fields 
in S7 + , Q\ and J) - 

p 7 ■ 

Let gl(x,y) correspond to y G 0^ and x G IX,-. Reciprocity implies gj(x.,y) = —gj(y,x). The 

integral representations (|49j) . (|5()|) and (j. r )l(l can finally be written in the condensed form (recalling 
that tt = ft + UfiiU fi -): 



p(y,«)= / 5o+(y,x) s °(x)^(x)+ / 5l (y )X )[(A ; (x) 2 -(fc 1 ) 2 )-^-v 



p(x)dn(x); Vy G ft 
(52) 



3 Comments on the integral representation of the field 



Eq. (|52|) can be written as: 



P s (y) :=p(y) - / 5o+(y,x) s °(x)^(x) 

5i(y,x) 



(Mx) 2 -(^) 2 )-^-V 



p(x)dfi(x); Vy G , (53) 
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wherein p s (y) is the field scattered by the inhomogeneity of the slab. 

This can be compared to the more common formulation employing the free-space Green's func- 
tion (G°(x,y)): 



p (y) 



G°(y,x)s°(x)dfi(x) 



+ 



/ G°(y,x) 



(kwy - (k»y - -z. ■ v 



p(x)dn(x); VyGO, (54) 



wherein p d (y,ui) is the field diffracted by the entire inhomogeneous slab (including the slab itself 
and its inhomogeneities) . The SGF formulation thus appears to be more suitable than the FSGF 
formulation, because the scattered field p s (y,u>) accounts at the outset for more of the physics of 
the interaction of the obstacle with the incident wave than p d (y). 

It can be shown that the neglected field is generally smaller in the SGF formulation than in the 
FSFG formulation. Effectively, when the SGF formulation is employed, the zeroth order Born ap- 

2 n.U2\ V P 



proximation consists in neglecting / gi(y,x 

/Hi 

,0 



(fe(x) 2 - (fc 1 ) 2 ) 



V 



P 



p(x)dO(x) compared to 



5 +(y ; x )s (x)dO(x) whereas when the FSFG formulation is employed, the zeroth-order Born 

s 

0+ 

approximation consists in neglecting J*^ G°(y,x) (fe(x)) 2 — (A; ) 2 — ^£ . V p(x)dO(x) in compar- 
ison to L G°(y,x)s°(x)dO(x). 

0+ 

The use of the SGF includes some multiple reflections, while the FSGF formulation combined 
with the Born approximation does not apply to high contrasts, because the employed linearization 
tacitly precludes multiple reflections. 

Finally, when the density is constant, the contrast component of the kernels of both formulations 
reduce to 



=(Z2) 



+ ia(w,x 2 ) 



kJ 



(55) 



wherein a(u,X2) is the absorption coefficient and j = for the FSGF and j = 1 for SGF. It has 
been shown in that the Born approximation is reasonable if the phase shift introduced by the 
inhomogeneous medium is less than it, i.e., weak and smooth heterogeneities of simple shape. The 
shift depends not only on the size, but also on the kernel, (eq. I55j) . i.e., on the frequency, on the 
absorption, and on the contrast between the two phase velocities. The use of the SGF, when the 
initial configuration is a homogeneous slab filled with a fluid-saturated porous material, allows us: 
i) to reduce the frequency dependence of the kernel, ii) to reduce the kernel itself by taking into 
account a phase- velocity that is closer to that of the host medium and also by taking into account 
the absorption (dissipation) of the material, and iii) to provide more accuracy, in the sense that 
on the one hand, the specific Green's function already accounts for dissipation and for some of the 
geometry of the problem and, on the other hand, the approximation of the field in the integral 
is more realistic than when the FSFG is employed. The usual way (i.e. when the FSFG is used) 
to avoid the problem induced by the absorption consists in adding some dissipation term in the 
approximated field in the integrand (Modified Born approximation), but not by acting directly on 
the kernel of the integral. Methods such as the distorted Born approximation, whose convergence 



12 



analysis has been carried out in [13], also allow to consider objects with larger contrast, but by acting 
only on the constrast function, i.e., without introducing additional effects on the approximated field 
in the integrand and on the Green's function used in the formulation. 

The SGF domain integral formulation thus allows the elimination of some of the disavantages 
of the FSGF domain integral formulation. This is obtained by acting on the kernel, the Green's 
function and the approximated pressure field, contrary to other methods employing the FSFG 
which act only on one or two of the components of the integrand. 

The combined effects of this action is to allow us to define and implement an iterative scheme, 
starting with the zeroth-order Born approximation and using the SGF formulation, to solve wave 
propagation problems involving a medium, whose components have high constrasts and in which 
there exist abrupt heterogeneities, so to consider objects with larger constrats, with respect to the 
surrounding medium, than would be possible with the conventional FSGF formulation. 

4 Specific ingredients of the computational procedure for the pre- 
diction of the field scattered by an inhomogeneous porous slab 
solicited by a plane incident wave. 

We now adapt the previous analysis to the determination of the field scattered by an inhomogeneous 
slab (the direction of inhomogeneity being X2) solicited by an incident plane wave. 

This type of incident wave is associated with s° = 0, so that it would appear that there is no 
solicitation in the above equations. Nevertheless, for an incident plane wave initially propagating 
in 0,q+, the integral over T^, l)39jl . does not vanish. It follows that the term corresponding to the 
solicitation takes the form of p 0+ (y), p 1 {y) and p° (y), which are the responses in the subdomains 
Oq+j ^1 an d f2 - respectively (i.e. the zeroth-order Born approximation) to an incident plane wave 
propagating initially in Q 0+ given by 



wherein we take (in the computations) uq = 100kHz to be the central frequency of the source 
spectrum. 

The zeroth-order Born approximation is given in appendix iBl 

4.1 Application of the first-order Born approximation in the SGF formulation 

We give here the explicit form of the first-order Born approximation within the framework of the 
SGF formulation. 

Remark: As a consequence of the separation of variables, all pressure fields can be written in 
the form 



(56) 



wherein k\ = k°s'm6 l , k^ % = k°cos9 l and 6 l the angle of incidence with respect to the +X2 axis. 
The spectrum of the incident takes the form of a Ricker-like wavelet of the form: 




(57) 




(58) 
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4.1.1 First-order Born approximation in the SGF formulation for y £ Qq + 

When y G fi +, 9i(y,x,^) = 5f? + (y,x, u), so that 

1 dp(x 2 ) d 



p(y,u) -p + (y,oj) 



/H-oo r-a 
/ 5? + (y,x) 
-oo J b 



{k(x 2 f - (k 



1\2\ 



• le lik iyi +ik°(y 2 -a)] a l ( a l cog ^1 ^ 2 _ b ^ 



ia° sin ^ h> 



p{x 2 ) dx 2 dx 2 
{k\ (x 2 -b)))dk 1 



p x (x) dxi dx 2 



2vr (2a°a 1 cos {k\l) - i ((a ) 2 + (a 1 ) 2 ) sin {k\l)) 



k 2 



{k{x 2 f - (k 



1\2\ 



1 dp(x 2 ) d 
p(x 2 ) dx 2 dx 2 



p 1 (x 2 ) 



+oo 



,i(fei-fei)^i dxi 



By making use of the identity (|12j) . ((ST?)) becomes 



(59) 



p(y,u) -p 0+ (y,uj) 

■jglifciw+i^-'dtt-o)]^,* ^ a i,< cos ( X2 _ 6 ^ _ ia o sin ^ 2 _ b ^j\ 



k 2 ' % ( 2a°> i a 1 > i cos ( k L /lj - i ((a -*) 2 + (a 1 -') 2 ) sin ( 
Introducing the expression of p 1 ^) from (|94|). and after expanding, we get: 

2ie^ fc i ?/1+ ^2' I (y2-2a)] Q ,l,j Q ,0,i 



p 1 (x 2 ) dx 2 . (60) 



( a l,i)2 _ ( a °.<) 2 . 



fc^ 1 ( 2a > i a 1 >* cos Mfc^lJ - i ((a -*) 2 + (a 1 -*) 2 ) sin ( 
( a l,i)2 + (a *) 2 x 



(A; 1 ) 2 ^ X (x 2 )dx 2 + V " 7 ^ V " ; (A: 1 ) 2 / X (x 2 ) cos ( 2fc^ (x 2 - b) ) dx 2 + 

itWlfe^ f i ^fr) cos f 2k 1 / (x 2 - b) ) dx 2 - 
Jb P{x 2 ) dx 2 1 - ' 



ia 1 -^ '^ 1 ) 2 y x(^2) sin [2k/ (x 2 - 6) J <*r 2 + 

(aW > 2 + ( «°" )2 r lags!*, (^» ( x 2 - 6) ) 

b p(x 2 ) <9x 2 V 



■ A, o 



(61) 



where x( x 2) = ( — tttto 1 I is the contrast function. 

V (* ) / 

We define the average value, cosine transform and sinus transform of a function f(x 2 ) = 
h(x 2 )U(b < x 2 < a) (wherein U(b < x 2 < a) is the so-called gate function and I = a — b) by: 
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< /(x 2 ) >=< h(x 2 ) >= J h(x 2 )^y- 

/oo 
f(x 2 )cos(qx 2 )dx 2 
-oo 
poo 

TF s (/(x 2 ),g) = / f(x 2 )sin(qx 2 )dx 2 



respectively. 

Eq. (|61j) can then be written in the form: 



2i e [i fc ^i+i fc 2 ,l (y2-2a)] Q ,i,i a o, l 



p(y,u) -p 0+ (y,w) 



(a 1 ' 1 ) 2 - (a '*) 2 



k 1 / ( 2a°> i a 1 > i cos ( ) - i ((a -') 2 + (a 1 .*) 2 ) sin ( k/l 



hi 



hi- 



(k ) I < X (x 2 ) > + 



■(fc A ) z TFc ( X(X2 ~ b),2k? ) + ia^a^TFc ( -^fo - 2# 



ia 1 '^ '*^ 1 ) 2 !^ (x(x 2 - 6), 2k£ t ) + 



i,A , (a M ) 2 + (a°-*) 2 , MrT ^ / 1 9p 



l.i 



^TFs(-^(* 2 -&),2£V 



(62) 



(63) 



The first-order Born approximation of the reflected field in the SGF formulation involves the av- 
erage value of x{ x 2) an d both the cosine and sinus transform of x( x 2~ b) and of —. — ^ 2 -, 

p(x 2 - b) dx 2 

while the first order Born approximation in the FSGF formulation involves only the Fourier trans- 
form of x( x 2 — b) and of — — ^ 2 defined by reference to the material parameter of the 



host |22|. 



p(x 2 - b) dx 2 



4.1.2 First-order Born approximation in the SGF formulation for y G fii. 

When y G Jli, <?i(y,x, u) = g\(y,x,u>), so that, proceding as previously, we get 



p(y,u) -p 1 (y,u) 



si(y. x ) 



p x (x) dxi dx 2 



/ 1 c te?/i+ifco''l?;2~a:2l1 1 jjklyi x 



a hif _ («o,<)2) cos fc M fa + X2 _ a _ & ) + e ifc 2 ' J 



a 



) z cos ffcj'* (y 2 - x 2 ) 



2k 1 / ( 2a°> i a 1 >< cos ( k/l ) - i ((a ^) 2 + (a 1 ' 4 ) 2 ) sin ( k/l 



hi 



{k(x 2 f - (k 



1 ^ 1 dp(x 2 ) d 

p{x 2 ) dx 2 dx 2 



p (x 2 ) dx 2 . (64) 
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4.1.3 First-order Born approximation in the SGF formulation for y £ &o-- 

When y G Q -, 5i(y,x,^) = g\ (y,x,u/), so that 



p(y,w) -p° (y,u) 



oo J b 



5?"(y,x) 



(/c(* 2 ) 2 - (k^) 



1 dp(x2) d 



p{x 2 ) 8X2 8X2 



p 1 (x) <ixi 



i e [ik\yi+ik 2 '\b-y 2 )] a i,i ( a i,i cos ^ ( -x 2 )) -ia^sinffc. 



c 2 ' (o-a? 2< 



( 2a°> i a 1 - i cos ( k l /l ) - i ((a -' ) 2 + (a^) 2 ) sin ( /c 2 M Z 



(fe(* 2 ) 2 - (fc 1 ) 2 ) 



1 dp(x 2 ) d 



hi 



p{x 2 ) dx 2 dx 2 



p 1 (x) d,X 2 ■ 



(65) 

Introducing the expression of p 1 (x) from (|94jl , and by making use of the definition Ijfi2|) , the 
previous equation can be written in the form: 



p(y,a/)-f> (y,w) 



1 i 



2a°' i a 1 ' i cos I fen' Z 



i((a°. i ) 2 + (a 1 - i ) 2 )sin (jfcj'\ 



(a 1 '*) 2 + (a -* ) 2 
2 

la^a^cos (fe^z) 
( a i,i)2 _ ( a 0,i)2 



cos ( fc^Z ) - ia u 'V' J cos ( fe 2 M Z ) ) l{k l f < X (x 2 ) > 



-1\2 



ft 



l,i\2 



cos ( k\ % l } ) fco''Z < 



cos \ k\ ,% l 

.,1,^2 



X ) 2 TF C ( *(x 2 - 6), 2^'* ) + sin ( fc^Z ) (fe i ) 2 TF s ( X (z 2 - b),2k 



hi 



l,i, 



1 d/?(x 2 ) 

p(x 2 ) cte 2 



> + 



l.i 



J\2r 



hi 



( a l,i)2 _ ( a 0,i)2 



cos ( jfeJ'V ) fc^TF s ( ~-g-(x 2 - 6), 2£; 2 



1 dp 



sin ( k\' % l ) fe 2 '*TFc 



- — (s 2 - 6),2A; 2 ' 
pox 2 



. (66) 

1 

This equation involves the average values, the cosine and sinus transform of both y(x 2 ) and — — — . 

p ox 2 

Compared with the formulae (|63j) . the transmitted field involves the additional term corresponding 

1 dp 

to the average value of — — — . 

pox 2 

4.2 The iterative scheme for solving the direct problem 

As pointed out previously, our aim is to define an iterative scheme to solve the direct problem of 
the diffraction of an incident plane wave by a heterogeneous porous slab. This would be of great 
interest for both the direct and inverse problems, due to the possible increased accuracy it can 
enable with respect to both the zeroth- and first-order Born approximations (in both the SGF and 
FSGF formulations). 

We want to compute the total fields in f2 + an d r2 — . These problems being formally similar, 
we will only detail the computation of p(y) ; y 6 Q + . 

Let p 0+ ^\y) and p 1 ^(y) designate the j-th iterates of the pressure fields in Q + and fii 
respectively. The iterative scheme proceeds as follows: 
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• Calculation of p 0+ ^(y) through (|6Uj). corresponding to the application of the Born approxi- 
mation in the SGF formulation. 

• Calculation of p 0+ ^\y) for j > 1. 

More specifically, we first have to calculate the pressure field in Qi by means of 



1(0), 



ie 1 **^ 1 ((a 1 -*) 2 - (a ' 4 ') 2 ) cos k\ % (y 2 + x 2 - a - b) + e 



2k 
i,t 



ifcjj/i+ifc 2 ' l |?/2-a;2| 



+ 



a 1 ' 4 ) cos ( (y 2 - x 2 ) 



{k{x 2 f - (fc 1 ) 2 ) 



2^ (2a°- i a 1 ' i cos (fcj'*z) - i ((a -') 2 + (a 1 '*) 2 ) sin 
1 dp(x 2 ) d 
p(x 2 ) dx 2 dx 2 



p 1{j - 1) (x 2 )dx 2 



wherein p 1 ^°^(y,a;) is the expression given in (j94j) . 

Once a new p^\y, oS) is evaluated, one computes a new p° ^(y) by means of the relation 



(67) 



p 0+ ^(y,a;)-p 0+ ( )( y , W )« 

"ieiM»i+i*2'*(jo-«) a i.* ( a 1 ' 4 cos ( Z^' 1 (x 2 - 6) ) - ia u sin ( /c^ (x 2 - b) 



k¥ ^aW* cos (kYl) - i ((a - 4 ) 2 + 



n 



V 7 p(x 2 ) <9x 2 9x 2 



^(xa)^. (68) 



Remark: Another scheme is the iterative calculation of p l ^\y,uj) and the subsequent computa- 
tion of the reflected field p 0+ ^(y, oj). 

The differentiation of p l ^\ which is a particular feature of our method, is carried out analyti- 
cally for j = by means of (|95j) . and numerically for j > 1 using the finite difference scheme: 



9 ~im (x s^ P m (i + l)-P m (i) 

dx 2 p [ 2) ~ x 2 (i + i)-x 2 (i) 



(69) 



a 



The computation of — — p 1 ^\x 2 

dx<2 X2 =X 2 (N)=a 

1 dp l< j\x 2 ) 

this derivative by using the fact that 



p(x 2 ) dx 



cannot be carried out in this manner. We approximate 
is conserved, so that 



_d_ 

dx< 



■P 1(j) (x 2 ) 



p( N ) 9 ,m 



X2=X2{N)=a p(N-l) dx 2 



p^\x 2 ) 



x 2 =X 2 (N-l) 



(70) 



wherein N is the number of discretisation points used to performed the calculation. 
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5 Outline of the numerical procedure 

We focus on the response of a double layer (each layer being homogeneous) porous slab (called layerl 
and layer2), considered to be a single inhomogeneous slab. We assume that the medium in the slab 
responds to a solicitation as does an equivalent fluid (i.e., this is the rigid- frame approximation). 

In an equivalent fluid medium, the appropriate conservation of momentum and constitutive 
relations take the form: 

u? V + 1 V • ( — L-Vp ) = (71) 

wherein 

p e (X, LO) = p e (x 2 ,Lo) = - L —z r 1 + 1 F{x 2 ,u) 

4>{x 2 ) V w 

1 1 7^0 (72) 



0(x 2 )( 7 -( 7 -l)(l + i^lG(x 2 ,PrM 

with w c {x 2 ) = - 2 ^f^ 2 ^ and G(:r 2 ,Pr 2 u;) 3 , F{x 2 ,u>)'2A\ being two relaxation functions given 
p f a{x 2 ) 

by 

cv ^ /i ■ 47 ? p / q 00 (x 2 ) 2 

F ^ -> - v 1 - -ggggggg^ (73) 

G(X2 ' Pr W) = V V(x 2 )^(x 2 ) 2 A-(x 2 ) 2 Pr W 
The chosen profile of porosity 4>{x 2 ), A(x 2 ), A'(x 2 ), a{x 2 ) and a{x 2 ) is presented tablc[5J 





4> 


Too 


A 


A' 




Thickness 






{(im) 


{(Jim) 


(A^.m" 4 ) 


{mm) 




Layer 1 


0.96 


1.07 


273 


672 


2843 


7.1 


Layer 2 


0.99 


1.001 


230 


250 


12000 


10.0 



Table 1: Properties of the two-layer medium studied. 

The inhomogeneous porous slab is included between b = — 10 x 10 _3 m and a = 7.1 x 10 _3 m. 
The contact surface between the two homogeneous porous sub-slabs, is located at x 2 = 0m. 
Special attention must be paid to: 

• the discretisation of x 2 in order to correctly model eventual jumps, 

• the modeling of the jump; as pointed out in appendix ^ the spatial dependence of the 
density p(x) can lead to meaningless integrals, especially when this parameter presents some 
discontinuities. To avoid this problem, we will consider such jumps to be well-approximated 
by the continuous function 

F(s 2 -e)«i(l + erf(£^)) (74) 
where e is the location of the jump, s the slope of the smooth jump and erf the error function, 
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• the determination of the parameters and /c 1 (a;) filling the initial homogeneous slab. 



5.1 Choice of the discretization step. 

Because of the necessary correct modeling of the 
at both the location of the step and the sides of 
increase of the point density at these locations, 
both sides of a jump. In our computations, A is 



continuous steps, making use of the formulae (|74|) 
the slab, we consider a logarithmic scale, with an 
This logarithmic scale occurs over a width A on 
chosen equal to 8 x 10 -6 . 



5.2 Choice of the modeling of the jumps 

To model the jumps at x 2 = 0, we define a function £ such that : 

C(*2) = Ci + ^y^(l + erf(^p)) (75) 

where ((x 2 ) can be 4>(x 2 ), A(x2), \'(x 2 ) : a(x 2 ), or a(x 2 ), and the indices 1 and 2 refer to the values 
of the parameter ( of the homogeneous layer 1 or 2. The quantities p(x 2 ) and k(x 2 ) are then 
computed. 

Once p 1 and k 1 are determined, in order to take into account the jumps at both (or at least 
one) sides of the entire slab, we compute: 

p{x 2 ) =p 1 + (p{x 2 ) - p l ) (erf(^^) - erf(^-p) - l) . (76) 

Thus, on both sides of the slab we model the "half" jump from p 1 to p(x 2 ) using a half of the erf 
function (compared to the jump inside the entire slab). This constitutes a better fit of the real 
jump. 

In all our computations, the parameter s is chosen equal to 2 x 10~ 6 . 



5.3 Choice of parameters p l (uj) and k l {u). 

The purpose of the SFG is to reduce the kernel of the integral ()53|) compared to the kernel of 
the integral Q54|) in the FSGF formulation. Because of the spatial dependence of the density, the 
integral (|53j) can be split into two integrals whose respective kernels are: 

(fe(x) 2 - (k 1 ) 2 ) andi|^. (77) 
p 0x2 

The easiest way to reduce these kernels would be (referred-to as choice 1.), all the characteristic 
parameters of the slab being known, to consider the average value of p(x 2 ,u>) and k(x 2 ,uj) over 
x 2 £ [b,a], as shown figureOJfor $t(p(x 2 ), 175kHz). 

Another choice (referred-to as choice 2.), which can be more convenient, consists in taking p l (oj) 
and A: 1 (w) equal to the minimal value of p(x 2 ,u) and k(x 2 ,u>) over x 2 G [b, a], as shown figure^] 
This choice would normally lead to the disappearance of the remanent density (i.e. equal to p 1 
whose values is larger than p(b)) at x 2 = b, figure El Another advantage of this choice is the 
reduction of the interval of integration , the first kernel vanishing over a part of this interval. 
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Figure 3: Real part of the density profile cor- 
responding to p 1 (175kHz) chosen as the av- 
erage value of p(x2, 175kHz) over X2 £ [b, a]. 



Figure 4: Real part of the density profile cor- 
responding to p 1 (175kHz) chosen as the min- 
imal value of p(x2, 175kHz) over X2 E [b, a]. 



To give an idea of the accuracy of the method, we introduce the following measure of the 
quadratic error, calculated for the j-th iteration: 



E 



Jo [ptmm( x ^) Psgim( x ^)) dt 
Jo (ptmm( x '*)) dt 



(78) 



wherein p°gQ IM (-x.,t) is the reflected pressure as computed by our Specific Green's Function based 
Iterative Scheme (SGIM), andp^M A// (x, t) the reflected pressure as computed by the classical Trans- 
fer Matrix Method (TMM), appendix The quadratic error corresponding to our computations 
is given figure 
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choice 2. 
choice 1 . 



7 8 9 1 



Iteration Iteration 
Figure 5: Evolution of the quadratic error as a function of the number of iterations. On the left: 

7T 

angle of incidence 0. On the right: — . 

3 

These experiments show that the correct choice of p 1 ^) and k (u) is indeed to consider the 
average value of p(x2,u) and k{x2,oS) over X2 £ [6, a]. This choice leads to a quicker and better 
convergence than the one obtained by the choice of p 1 ^) and /c 1 (cj) as the minimum of p(x2,u>) 
and k(x2,uj) over X2 £ [b,a\. 

For both choices of these parameters, after a certain number of iterations, the SGIM results are 
the same as the classical TMM results, as shown figure El for example, when choice 1 is made. 




t (|1 S) t (|1S) 

Figure 6: Reflected pressure as computed by the classical Transfer Matrix Method (TMM) -dashed 
curve- and as computed by our Specific Green's Function based Iterative Method (SGIM) when 

7T 

choice 1. is made. The angle of incidence on the left and — on the right. 

3 

1 dp 

Remark: For both choices, < — — — >, involved in the calculation of the first order Born 

P OX2 
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approximation of both the reflected and transmitted fields, vanishes, i.e. 

ra 1 dp 



1 dp 

< 'IT- >z 
p 0x2 



-dx<} 



lb 9 dx 2 

because p(x 2 ) = p l {oj) for both x 2 = a and x 2 = b 



a dln(p(x 2 )) d ^ = ^ /p(a) 
dx 2 2 \p(b) 



(79) 



Remark: Other choices are possible, such as the one leading to the disappearance of the av- 
erages < x( x 2) > involved in the calculation of the first order Born approximation of both the 
reflected and transmitted field. Consider < x( x 2) >■ 



< X{ 



( w f7(W ,\ dx 2 1 r 2 dx 2 

iX2) >= I Iw - 1 )— = w?k {k{X2)) — 

which vanishes only if k l = \ 



1 



(80) 



, dx<2 



(k(x 2 )) 2 ——. This choice corresponds to the particular case in 
lb l 

which k(x 2 ) is such that the Schwartz inequality is satisfied: 



(81) 



J^) </w^. 

6 Results and discussion 

We first present results, as calculated by the iterative scheme (FGIM) initialized with the zeroth- 
order Born approximation arising from the integral formulation incorporating the free-space Green's 
function, to emphasize the fact that this method does not converge in all cases. 
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Figure 7: Reflected pressure as computed by the classical Transfert Matrix Method (TMM) and as 
computed by the classical free-space Green's function based iterative scheme (FGIM). On the left : 

7T 

the incidence angle is ; in the middle : the incidence angle is — and on the right : the incidence 

6 

angle is — . 

For small angles of incidence, the usual FGIM converges, figure |H] and [7J but slower and with 

7T 

less accuracy than our SGIM (figure EJ. For large angles of incidence (in our example, — ), the 
usual FGIM strongly diverges, figure [71 while our method still rapidly converges. This is probably 
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caused by the fact that the first iteration is far from the exact solution (figure |HJ). The translation 
of this divergence can be appreciated in (figure 



Incidence angle 

Inicidence angle rc/6 

Incidence angle Jt/3 











1 2 3 45678 910 15 20 25 

Iteration 



Figure 8: Evolution of the quadratic errors as a function of the number of iteration for various 
angles of incidence in the FGIM method. 

All the following computations are carried out with characteristic parameters p 1 (w) and A: 1 (a;) 
filling the initially-homogeneous slab chosen as the average, over a < x% < b, of p(x2,uj) and 
k 1 (x2,to) respectively. 

We define a convergence criterion via the quadratic difference between two iterations i and j : 

A * = — 5 • (82) 

We found empirically that a convergence criterion C(i) of the form 

C(i) is true when ^2±M. < j x iq- 6 (83) 

£>2,1 

gives good results as shown figure El and E3 Our method yields solutions that are close to those 
of the usual TMM both in transmission and reflection in both the frequency and the time domain 
for several angles of incidence. This validates the method employing the SGF, combined with an 
iterative scheme initialized by a zeroth-order Born approximation. 

The time history of the reflected pressure figure is of particular interest for the demonstration 
of the accuracy of our method. We can clearly distinguish, for both angles of incidence and 

7T 

— , the three reflections of the incident wave on the three interfaces of our canonical configura- 
3 

tion. The zeroth-order Born approximation, i.e. corresponding to the homogeneous fluid-saturated 
porous slab, formally accounts for two of them. For the first and third reflection, the amplitude 
matches correctly with our convergence criterion. The second reflection, which comes from the 
inhomogeneous slab, matches in both amplitude and time of arrival. 

The time history of the transmitted pressure figure E3 contains only one peak due to the fact 
that absorption within the slab attenuates the transmitted waves resulting from multiple reflection. 
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Figure 9: Reflected pressure. At the top: Convergence criterion; in the middle: the reflected 
pressure; at the bottom: the spectrum of the reflected pressure. On the left: incidence angle = 0. 
On the right: incidence angle 
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Figure 10: Transmitted pressure. At the top: Convergence criterion; in the middle: the transmitted 
pressure; at the bottom: the spectrum of the transmitted pressure. On the left: incidence angle 
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0. On the right: incidence angle = — . 



7 Conclusion 

A method, making use, in the domain integral formulation, of the specific Green's function (SGF), 
i.e., the Green's function of a canonical problem close to the original problem, for the resolution of 
problems of acoustic wave propagation in an inhomogeneous fluid medium (with spatially-varying 
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density and compressibility) was studied and implemented for the canonical example of plane wave 
solicitation of a double layer fluid-saturated porous slab (considered as a single inhomogeneous 
slab) in the rigid frame (equivalent fluid) approximation. 

A particular feature of our study is that we account for spatially-varying density, contrary to 
many authors who consider it to be constant. We also address the issue of the spatial differentiation 
of the pressure field at the boundaries of the inhomogeneity, which is carried by a finite-difference 
scheme for higher-than-zeroth-order Born approximations. 

Our specific Green's function iterative scheme, which is initialized by a zeroth-order Born ap- 
proximation, was shown to converge, contrary to the iterative scheme relying on the free-space 
Green's function, which is often found to be divergent. This improvement is due to the com- 
bined effects of the use of the SGF and to a better Born-like approximation of the field inside the 
heterogeneity. 

In our numerical examples, our method was found to converge within 5-7 iterations to the 
reference solution (obtained rigorously by the transfer matrix method), even for an abrupt hetero- 
geneity, and for various choices of the acoustic parameters filling the homogeneous slab supposed 
to be the initial configuration (canonical problem) for the SGF. 

The robustness of the method was also demonstrated. 

Our method thus appears useful for the resolution of inverse problems. In such a context, some 
information about the geometry and/or the mechanical properties of the objects one is looking for, 
is often known. The SGF is the device by which this information can be incorporated into the 
inversion procedure in a rational manner. 



A The wave equation in an inhomogeneous fluid medium 

A.l Solution of the direct problem involving both the pressure and its partial 
derivative 

Wave propagation, relative to an acoustic wave in an inhomogeneous fluid occupying a domain 
= $7o U ill (the homogeneous host medium occupies the domain while the inhomogeneity 
occupies the domain fii), is described by: 

V ' Vp+ ^-^^ Vp = / ° (x)s(x) ;VXG0 (84) 



wherein: c(x) = y - j — is the spatially-varying velocity, and k(x) and p(x) the spatially- 
varying compressibility and density respectively of the fluid. 

Applying the domain integral formulation with the usual free-space Green's function, leads to 
the domain integral representation of the total field 



p(y)=P*(y)+ / G°(y,x) 



UJ 



2 



c(x)S 



(k 



0\2 



P(y) - • Vp(x) ) dfi(x) ; Vy £ $1 (85) 



wherein p l (x) is the incident field. To obtain the field at an arbitrary point of space, p and Vp 
within Q\ have to be determined. This can be done by solving the coupled system of integral 
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equations: 

p(y)=p*(y)+ / G°( y ,x) 



c x 



0\2 



p(y)-^7^-Vp(x)] ( /0(x) : VyO, 



V y Ky) = V y p J (y)+ / V y G u (y,x) 
./Ql 



c x 



(A; 1 



P(x) 

" p(y) - • v p( x ) ) d ^( x ) ; v y^i 



P(x) 



(86) 

A vast literature exists on the subject of the numerical resolution of systems of domain integral 
equations [2T1 f27H IH1 fTTH IFH UK] . 

A. 2 Solving the direct problem via a single integral equation 

Another, perhaps simpler (although unsuitable in the inverse problem context) way, to solve the 
previous problem is to make the substitution 

p(x)=g(x)Vp(x) (87) 

whereby the following governing equation is obtained 



V 2 g(x) + 



uJ 



1 V 2 p(x) 3 Vp(x) • Vp(x) 



g(x) = pa (x)s(x) ; x G ft 



_c(x)2 2 p(x) 4 (p( x )) 2 
Using the free-space Green's function in the domain integral formulation yields the representation 

q(y) = q i (y)+ 



[ G°(y,x) 



, ,2 , lV-Vp(x) 3Vp(x)-Vp(x) 



q(x)dv(x) ; Vy G fi, 



.(c(x)) 2 VW ' 2 p(x) 4 (p( x )) 
from which is extracted a single integral equation for g(x) ; x G fii. 

Remark: The integral formulations ()86|) and (|89() are identical when the density is constant. 
A. 3 A canonical problem involving a density discontinuity. 

Let us consider the simple ID problem, depicted in figure II 1| of a plane wave striking a planar 
interface T, located at X2 = a, between two homogeneous media and fli. The normally-incident 
plane wave travels initially in f^o- The heterogeneity is supposed to be the domain Q±. 

In practice, this problem can be treated rigorously by the TMM method. However, when one 
attempts to solve it by the integral method, the medium filling fl\ must be dissipiative. 

Let us suppose that the pressure field p 1 in fii is known (for example, calculated by the TMM 
method). 
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a 

Figure 11: Configuration of a pianar interface between two domains. 

We introduce p(x 2 ) = p 1 + (p° — p 1 ) H{x 2 ) and k(x) = k 1 + (A; — A; 1 ) H(x 2 ), where H{x2) is 
the Heaviside function, and k J , j = 0, 1 the wavenumber in the domain Qj. 
Eq. @2J) splits into: 



p{V2)=p i {V2)+ I G (y 2 ,x 2 ){k 1 ) 2 p(y)dx 2 -{p°-p 1 ) 



u f a S(x 2 - a) dp(x 2 ) 



p(x 2 ) dx 2 



dx 2 ,\/yn (90) 



wherein 5 is the Dirac delta distribution. All the integrals involved in (|9U|) can be solved analytically 

1 dp(x 2 } 

(p(x 2 ) begin known by hypothesis). In particular, the function — — - — is continuous (i.e. the 

p{x 2 ) ox 2 

function is C°) at the interface T so that the second integral does not present any difficulties. 

The formulation involving the evaluation of p and of Vp allows us to take into account density 
discontinuities. Eq. lf%9"|) splits into 



g(y 2 ) = <f(y 2 ) + / G° (y 2 ,x 2 ) {k 1 ) 2 g (x 2 )dx 2 + 

[P°-P l ) f a S'(x 2 - a) 
2 i-oo p(.x 2 ) 
,0 „1\ 2 f -a 



q(x 2 )dx 2 - 



S^-P 1 ) f a 5(x 2 - a)5{x 2 - a) 

(P(x)) 2 



<{x 2 )dx 2 ■ Vy 2 G n • (91) 



The calculation of first integral presents no particular difficulties. Let us consider the second integral 
' "(x 2 )dx 2 , wherein S'(x 2 — a) is the derivative of the Dirac delta distribution. The 



-oo P( x z) 
use of the formula 



f(x 2 )5'(x 2 - a)dx 2 = - ^—*- 6(x 2 - a) 

ox 2 



(92) 



requires the function f(x 2 ) to be C 1 at x 2 = a, while the function — 2 | = ^ 2 \ is not continuous 

P(X2) p(x 2 )2 

at the interface T. Thus, this second term cannot be handled analytically. 
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Finally, consider the third term f - 2 - ^ -q(x2)dx2- The integrand involves the 

J-oo (p(x)) 

scalar product of two Dirac delta distributions 5{x2 — a)5{x2 — a) which is not defined |31H15lHl) . 
A Numerical approximation of this quantity exists, but otherwise it is meaningless 114. . 

The resolution of problems via a single equation is also of no practical use when the problem 
one is faced with involves density discontinuities. 

B Pressure field in the case of a macroscopically-homogeneous 
porous slab (zeroth-order Born approximation). 

By referring to ^Hj 5 one finds, that for plane wave solicitation in Qq, the pressure fields in Qq , f^i 
and Qq are: 



p (x, uj) = A l (u>)exp yik\x\ — ik 2 ' x 2 j + 

iexp (ik\x x + ik°/(x 2 - 2a)) sin (hfrH) ({a 1 ' 1 ) 2 - (c^) 2 ) 



A\u) ^- ~ T —^ J - v " 7 . , — = exp (ik\x x ) p° (x 2 ,lo) (93) 

2a 1 > i a°>* cos (k l 2 ' l Vj - i ((a 1 ^) 2 + (a ^) 2 ) sin (k^l) 



2 exp (ik\xi — ik^a 


a '* 


a 1 ' 1 


cos (k\ ,l {x2 — b) 


— ia ' 4 sin 


[ kl'\x 2 - &)) 




2a 1 > i a°> i cos (k 1 /^ 


-i((a 1 ' i ) 2 + (a c 


> i ) 2 )sin (k\> 


'0 





exp (i/sfa;i)p 1 (a;2,a;) (94) 



(9p 1 (x cj) 2 exp ^i/c|xi — ik^aj cP ,l k\' 1 —a 1 ' 1 sin ( k 2 ~' i (x 2 — b)\ — ia 0,i cos I /^(^ — b) 



• ' ' ' ~ ' 2 a M a o,i cos (fc^z) - i ((a 1 -*) 2 + (a ^) 2 ) sin 



2 exp ( i^xi - ifejj'* (a? 2 + ) a 1 '*" '* 



p° (x,u) = A\u) - 7TV77{ 77— „ 7~T7A = exp (ifc^i) p° (x 2 ,w) . 

(96) 



2a 1 ' i a°> i cos ( fc^J - i ((a 1 '*) 2 + (a 0,i ) 2 ) sin ( fc 2 M Z 
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C Pressure field in the case of a double layer macroscopically- 
homogeneous porous slabs. 

We use a separation of variables technique to obtain the field representations: 

p 0+ = A^f^e^** 1- **'*^ _j_ B 0+ e \[k\x 1 +k 2 ' l (x2-a)] 
pi _ gifcl^i ^l e -ifc2 ,l (^2-a) _|_ j^l e ikl'\x 2 -a) S ' 
p 2 = gifcjxi ^ A 2 e -lk¥(x 2 -b) + B 2 e lk 2 2 '\x 2 - 
p°~ = J^- e \[k\x 1 -kl'\x 2 -b)] 



(97) 



After introducing the fields expressions into the boundary conditions (continuity of the pressure 
and of the normal velocity), we multiply these relations by exp (— \K\Xi) and then integrate form 
— oo to +00 to obtain the matrix equation (solved numerically to get B 0+ and ^4° ) 



1 


-1 


-1 













a 1 '* 


-a 1 ' 4 
















































1 


1 


-1 











-a 2 ' 4 


a 2 '* 





/B 0+ \ 
A 1 
B 1 
A 2 
B 2 

\ A ' r J 



( -A i {u)e-' lk ^ a \ 
a 'M i (w)e- ifc °' v 







(98) 

h ■-) 



wherein I 1 and I 2 are the thickness of the layer 1 and the layer 2 respectively (tableEI) and a J ' 1 = —. 

p] 

j = 0+, 1,2,0-. 
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